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Certain numerical frameworks used for the evolution of binary black holes make use of a gamma 
driver, which includes a damping factor. Such simulations typically use a constant value for damping. 
However, it has been found that very specific values of the damping factor are needed for the 
calculation of unequal mass binaries. We examine carefully the role this damping plays, and provide 
two explicit, non-constant forms for the damping to be used with mass-ratios further from one. Our 
analysis of the resultant waveforms compares well against the constant damping case. 

PACS numbers: 04.25.D-, 04.25. dg, 04.25.Nx 



I. INTRODUCTION 

The ability to simulate the final inspiral, merger, and 
ring-down of black hole binaries with numerical relativ- 
ity plays a key role in understanding a source of 
gravitational waves that may one day be observed with 
gravitational wave detectors. While initial simulations 
focused on binaries of equal-mass, zero spin, and quasi- 
circular inspirals, there currently is a large efFort to ex- 
plore the parameter space of binaries, e.g. 043- ^ ^^y 
part of studying the parameter space is to simulate bina- 
ries with intermediate mass-ratios. 

To date, the mass ratio furthest from equal masses that 
has been numerically simulated is 10:1 These sim- 

ulations use the Baumgarte-Shapiro-Shibata-Nakamura 
(BSSN) formulation |lflilj| with 1-t-log slicing, and the 
f driver condition for the shift [l^ [l^l- In it was 
noted that the stability of the simulation is sensitive to 
the damping factor, 77, used in the F driver condition. 



(1) 



Here, /3* is the shift vector describing how the coordinates 
move inside the spatial slices, do = dt — f3^di, and I" is 
the contraction of the ChristofFel symbol, r* j,, with the 

conformal metric, 7-''"'. 

The standard choice for rj is to set it to a constant 
value, which works well even for the most demanding 
simulations as long as the mass ratio is sufficiently close 
to unity. In binary simulations, a typical choice is a con- 
stant value of about 2/M, with M the total mass of the 
system. This choice, however, leads to instabilities for 
the mass ratio 10:1 simulation Q, although stability was 
obtained for rj = 1.375/M. The value of r/ is chosen to 
damp an outgoing change in the shift while still yield- 
ing stable evolutions. As we will show, if r/ is too small, 
there are unwanted oscillations, and values that are too 
large lead to instabilities. By itself, this observation is 
not new, see e.g. [TB - flsj . The key issue for unequal 
masses is that, as evident from ([T]), the damping factor r] 
has units of inverse mass, 1/M. Therefore, the interval 
of suitable values for 77 depends on the mass of the black 



holes. For unequal masses, a constant r] cannot equally 
well accommodate both black holes. A constant damping 
parameter implies that the effective damping near each 
black hole is asymmetric since the damping parameter 
has dimensions 1/M. For large mass ratios, this asym- 
metry in the grid can be large enough to lead to a failure 
of the simulations because the damping may become too 
large or too small for one of the black holes. To cure this 
problem, we need a position-dependent damping param- 
eter that adapts to the local mass. In particular, we want 
it to vary such that, in the vicinity of the i^^ puncture 
with mass Mi, its value approaches I /Mi. 

A position-dependent 77 was already considered when 
the f driver condition was introduced [sl. [l9l - l22| . but such 
constructions were not pursued further because for mod- 
erate mass ratios a constant rj works well. Recently, we 
revived the idea of a non-constant rj for moving punc- 
ture evolutions in order to remove the limitations of a 
constant 77 for large mass ratios. In |23j . we constructed 



a position-dependent 77 using the the conformal factor, 
-0, which carries information both about the location of 
the black holes, and about the local puncture mass. The 
form of T] was chosen to have proper fall-off rates both 
at the punctures and at large distance from the binary. 
In , this approach was used successfully for mass ratio 
10:1. (We note in passing that damping is useful in other 
gauges as well, e.g. in |24| the modified harmonic gauge 
condition includes position-dependent damping by use of 
the lapse function.) 

In the present work, we examine one potential short- 
coming of the choice of H^, which leads us to suggest 
an alternative type of position-dependent r]. Using [23[, 
we find large fluctuations in the values of that rj, and 
this might lead to instabilities in the simulation of larger 
mass-ratio binary black holes. To address this, we have 
tested two new explicit formulas for the damping factor 
designed to have predictable behavior throughout the do- 
main of computation. We find the new formulas to pro- 
duce only small changes in the waveforms that diminish 
with resolution, and there is a great deal of freedom in the 
implementation. Independently of our discussion here, in 
[16| the stability issues for large rj are explained, and a 
non-constant rj is suggested (although not yet explored 
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in actual simulations), that, in its explicit coordinate de- 
pendence, is similar to one of our suggestions. 

The paper is organized as follows. We first describe the 
reasons for the damping factor and some of the reasons 
for limiting its value in SecUll In Sec. lIIIi we discuss some 
previous forms of 77 that have been used. We also present 
two new definitions and why we investigated them. In 
Sec. llVi we find that these new definitions agree well 
with the use of constant j] in the extracted gravitational 
waves for mass ratios up to 4:1. Finally, in Sec. [Vj we 
discuss further implications of this work. 



II. MOTIVATION 

In order to define a position-dependent form for rj, it 
is important to determine what this damping parameter 
accomplishes in numerical simulations. For this reason, 
we examine the effects of running different simulations 
while varying rj between runs. First we use evolutions of 
single non-spinning black-holes to identify the key phys- 
ical changes. Then we examine equal-mass binaries to 
determine specific values desired in 77 at both large and 
small radial coordinates. 



A. Numerics 

For all the work in this paper, we have used the BAM 
computer code described in [OjHEHBl- It uses the BSSN 
formalism with l-|-log slicing and F driver condition in 
the moving puncture framework [J [27| . Puncture ini- 
tial data [2^ with Bowen-York extrinsic curvature (29| 
have been used throughout this work, solving the Hamil- 
tonian constraint with the spectral solver described in 
[soj . For binaries, parameters were chosen using [3l| to 
obtain quasi-circular orbits, while the parameters for sin- 
gle black holes were chosen directly. We extract waves 
via the Newman- Penrose scalar ^'4. The wave extraction 
procedure is described in detail in [l^. We perform a 
mode decomposition using spin-weighted spherical har- 
monics with spin weight —2, as basis functions and 
calculate the scalar product 



sin6ld6'd(/3y;-^^'4. (2) 



We further split ^'4" into mode amplitude Ai„i and phase 
(j>im in order to cleanly separate effects in these compo- 



nents, Tex ■ '^4 = AimB 



In this paper, we focus on 



one of the most dominant modes, the I = m = 2 mode, 
and report results for this mode unless stated otherwise. 
The extraction radius used here is = 90 M. 



B. Single, non-spinning puncture with constant 
damping 

The damping factor, 77, in Eq. ([l}, is included to reduce 
dynamics in the gauge during the evolution. To examine 
the problem brought up in the introduction, we compare 
results of a single, non-spinning puncture with mass M. 
We use a Courant factor of 0.5 and 9 refinement levels 
centered around the puncture. The resolution on the 
finest grid is 0.025 M, and the outer boundary is situated 
at 256 M. Varying the damping constant between 0.0/M 
and 4.5/M, two main observations can be made. First, as 
designed, a non-zero rj attenuates emerging gauge waves 
efficiently. Second, an instability develops for values of ry 
that are too large. 

Figs, m and [H illustrate the first observation. Both 
figures show the x-component of the shift along the x- 
coordinatc using r] e {0.0/Af, 1.5/M, 3.5/M}. Apart 
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FIG. 1: The x-component of the shift, /3^, for a single non- 
spinning puncture of mass M at time t = 15.2 A'l. The three 
lines were taken for different values of the damping factor 77. 
The solid line (black) is for rj = 0.0/M. The dashed line (red) 
is for Tj — 1.5 /M and the dotted-dashed line (green) is for 
r] = 3.5/Af. This shows the beginning of a pulse in for 
smaller values of rj. 

from the usual shift profile. Fig. [1] shows the beginnings 
of a pulse in the 77 = 0.0/M case (solid line) at x w lOM 
after 15.2 M of evolution. Examining Fig. [51 where we 
zoom in at a later time, t = 30.4 M, one can see that 
the pulse has started to travel further out (solid line). 
Looking carefully, one can also see a much smaller pulse 
in the 77 = 1.5/M line (dashed). Lastly, by examination, 
one can find almost no traveling pulse in the 77 = 3.5/M 
curve (dotted-dashed line). The observed pulse in the 
shift travels to regions far away from the black hole and 
effects the gauge of distant observers. This might have 
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FIG. 2: The a:-component of the shift, (3^, for a single non- 
spinning puncture of mass M at time t = 30.4 M. The three 
hnes were taken for different values of the damping factor rj 
with the same line type and color scheme as in Fig. [1] Here 
it is clear a pulse radiates outward in the shift with smaller 
values of r/. 



undesirable implications for the value of such numerical 
data when trying to understand astrophysical sources. 

For values of 77 larger than 3.5/Af, an instability 
arises in the shift at larger radius. Fig. [3] shows the x- 
component of the shift vector using damping constants 
T] = 3.5/M (soHd line), -q = 4.0/M (dashed Hue) and 
77 = 4.5/M (dotted-dashed line). The plots show an in- 
stability in simulations with rj > 3.5/M developing in 
/3% which eventually leads to a failure of the simulations. 
Contrary to this, the simulation using rj ~ 3.5/M does 
not show this shift related instability. In test runs we 
found that by decreasing the Courant factor used, we 
could increase the value of the damping factor and still 
get stable evolutions. This agrees with where it was 
shown that the gamma driver possesses the stiff prop- 
erty, which limits the size of the time-step in numerical 
integration based on the value of the damping. 

Figures [1] [2 and [3] make clear how the choice of the 
damping factor affects the behavior of the simulations. 
The value we choose for 77 should be non-zero and not 
larger than 3.5/M to allow for effective damping and sta- 
ble simulations. The exact cutoff value between stable 
and unstable simulations is not relevant here since the 
position dependent form we develop in Sec. IIIII gives us 
the flexibility we need to obtain stable simulations. 
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FIG. 3: The s-component of the shift vector in the a::-direction 
for a single non-spinning puncture of mass M at times t = 
16.0 M and t = 30.4 M. The three different lines mark three 
values of the damping constant rj. The solid line (black) is 
for 77 = 3.5/M, the dashed line (red) for rj = 4.0/M and the 
dotted-dashed line (green) for 77 = 4.5/M. At t = 16 M, the 
simulation using 77 = 4.5/M develops an instability in the 
shift vector and fails soon afterward, the same happens for 
77 — 4.0/A/ at t = 30.4 Af. In the simulation using rj = 3.5, 
no such instability develops (not shown). 



C. Equal mass binary with constant damping 

To examine the effect of 77 on the extraction of gravi- 
tational waves, we compare the results from simulations 
of an equal mass binary with total mass M in quasi- 
circular orbits with initial separation D = 10 A/, using 
77 G {0.0/A/,0.5/M,2.0/M}. Again, the Courant factor 
is chosen to be 0.5 and we use, in the terminology of [iTj . 
the grid configuration x[6 x 56 : 5 x 112 : 6] with a finest 
resolution of 0.013 M. Here, the extraction radius is 
chosen to be 90 A/. 

For vanishing 77, we find a lot of noise in the the real 
part of the 22-mode of shown in the solid curve of 

Fig. [31 A small, but non-vanishing 77 suffices to suppress 
this noise, as seen in the dashed curve of this figure. The 
dotted-dashed curve in this plot is the result for using the 
value 77 = 2.0/Af. We sec a difference in time between 
peak amplitudes of the three curves due to the change 
of coordinates that the alternation of 77 introduces. We 
did, however, find that by decreasing the Courant factor 
those differences between peak amplitudes summarily de- 
creased. 

To understand the noise in the waves for 77 = 0.0/Af , 
we look at the shift vector at different times. The first 
panel of Fig. O shows the x-component of the shift over 
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FIG. 4: Real part of the 22-mode of over time for equal 
mass simulations using different values for t). The inset shows 
the full waveform until ringdown. The solid curves (black) 
are for 77 = 0, the dashed curve (red) mark r] = 0.5/Af and 
the dotted-dashed curves (green) are for r) = 2.Q/M. Without 
damping in the shift, the extracted waves are noisy at times 
when the amplitude is still small (black, solid curve). 



X, again for 77 e {0.0 /A/, 0.5 /A/, 2.0 /A/}, shortly after the 
beginning of the simulation. The fourth panel shows the 
same at a time shortly before the merger, and the two 
panels in the middle represent intermediate times. We 
see clear gauge pulses in the earliest time panel for all 
three curves. We also observe the amplitude of this pulse 
decreasing with increasing j]. As time goes on, the gauge 
pulse travels outwards as in the case for a single puncture 
in section IIIBI For vanishing 77 (solid line), the shift 
becomes more and more distorted, and the distortions 
do not leave the grid. For non-zero 77, the amplitude of 
the gauge pulse decreases when traveling outwards, and 
the shape of is not distorted. There is, compared to 
77 = 2.0/M (dotted-dashed line), only a small bump left 
in the 77 = 0.5/A/ case (dashed line), that changes its 
shape slightly during the simulation, but does not travel 
to large distances from the punctures. The coordinates 
are disturbed in the case where no damping is used, and 
thus the noise in rex-Re{^'4^} is not surprising. 

In this series, using a Courant factor of 0.5, we only 
obtained stable evolutions for 77 < 3.5/Af which agrees 
with the limits found in section lTl Bl If we chose the value 
of 1] too large, the same kind of instability in the shift 
vector wc found there develops in the equal mass case and 
the simulations fail. The failure occurs relatively early, 
before 50 M of evolution time, whereas the stable runs 
lasted about 1200 Af , including merger and ringdown (we 
stopped the runs after ringdown). 



III. POSITION-DEPENDENT FORMS OF 77 



In section lTl Bl we saw that a sufficient level of damping 
is needed to limit gauge dynamics, and too much damp- 
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FIG. 5: a;-component of the shift vector, for three differ- 
ent choices of r\ at four different times during the simulation. 
The physical system is the same as in Fig. [l] The merger 
takes place at approximately t — 1000 Af. In the 77 = 0/M 
case (black, solid curve), the shift vector is not damped and 
therefore, a pulse travels outwards and distorts the shift over 
the whole grid. The amplitude of this pulse is considerably 
damped when using a non-vanishing rj and therefore the dis- 
tortions are reduced. For 77 — 0.5/A// (red, dashed curve), 
there are still small bumps traveling out which are reduced 
by using 77 = 2/Af (green, dotted-dashed curve). 



ing can lead to numerical instabilities. In section III CI 
we saw the positive effect that sufficient damping has on 
the resultant waveform for equal mass binaries. While we 
still need damping in the gamma driver in the unequal 
mass case, a constant value may not fulfill the require- 
ments of limiting gauge dynamics and permitting stable 
evolutions. Rather, we need a definition for the damping 
that adjusts the value to the local mass-scale. 

We will examine definitions, that naturally track the 
position, and mass of the individual black holes. The 
choice of rj should provide a reasonable value both near 
the individual black holes, and at large distance from 
the binary. We will start by examining some previous 
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work, that has used non-constant forms of the damping 
parameter, and why it may be necessary to use other 
formulas. We will then present the two new formulas for 
T], which we designed for this work. 



A. Previous dynamic damping parameters 




800 



A position dependent damping was introduced some 
years ago by the authors of |20|, and was later used in 
[2l| . That formula reads 



^ — '?punc 



^punc ^cjo 



being constants, and assuming ip 



(3) 



vth 



with ?7pu„c, ?7oo 

1 + Mi/(2ri) + M2/(2r2) (r^ is the distance to the i 
puncture) . This formula was used to damp gauge dynam- 
ics while using excision for equal-mass head-on collisions. 
It has since been found that using the moving puncture 
framework allows for constant damping in the approxi- 
mately equal mass case. Wc are looking for a formula 
which is suitable for the quasicircular inspiral of inter- 
mediate mass-ratio binaries. 

Previously [l^, we used the formula 



r]{r) = Ro 



(l-^/'-2)2 



(4) 



for determining a position dependent damping coeffi- 
cient instead of using a constant rj. With i?o taken to be 
a unitless constant, it can be seen that Eq. ^ has units 
of inverse mass. The dependence on the BSSN variable, 
Tp, naturally tracks the position, and mass of the black 
holes. The application of Eq. (j4]) gave good values for 
the damping both at the punctures, and at the outer 
boundary, and was even found to somewhat decrease the 
grid-size of the larger black hole. The latter point could 
have positive effects on how the individual black holes 
are resolved on the numerical grid. It even had the addi- 
tional effect of keeping the horizon shapes roughly circu- 
lar, even close to merger - something that doesn't hold in 
the constant rj case. Most importantly, the simulations 
remained stable, without significantly changing the grav- 
itational waves. The formula was later used successfully 
for the 10:1 mass-ratio in 

Despite all this, Eq. (U) provides reason for concern. 
Fig.[6]shows the form of 77 using Eq. (|4]) for a non-spinning 
binary of equal mass in quasicircular orbits starting at a 
separation of I? = 10 M at four different times in the 
simulation. As can be seen, noise travels out from the 
origin as time progresses. This leaves steady features on 
the form of 77 which could spike to higher and lower val- 
ues than the range determined in Sec. Ill Bl Additionally, 
these sharp features may lead to unpredictable coordi- 
nate drifts, and could, in some cases, affect the long-term 
stability of the simulation. 

To illuminate the origin of the disturbances in 77 (r), 
we looked at the development of ri{r) in simulations of 
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FIG. 6: Damping factor, 77, along the a;-axis using Eq. Q. 
The simulated configuration is an equal mass binary with ini- 
tial separation Z) = 10 M and orbits lying in the [x, j/)-plane. 
Shown are four different times during the simulation. 



a single, non-spinning puncture, and a single, spinning 
puncture (S'z/M^ = 0.25). The result for the spinning 
case is plotted in Fig. [7] at two different times over the 
X-axis. Again, we see a pulse traveling outwards. Only 
this time, it does not leave much noise on the grid. The 
fact that this pulse travels at a speed which is roughly 
1.39 (in our geometric units where c = G = 1) in both 
the spinning and non-spinning scenario indicates that it 
is related to the gauge modes traveling at speed in 
the asymptotic region where a ~ 1 (see |32| and [Toj for a 
discussion of gauge speeds). In contrast to gauge pulses 
in the lapse, a, or shift vector, /?', the pulse in rjix) is 
amplified as it walks out. We found the same result in 
the single puncture simulation without spin. We believe 
the reason for this behavior is that as the distance to the 
puncture increases, the conformal factor, 7/), gets closer to 
unity. Therefore, the denominator in Eq. Q approaches 
zero, and the gauge disturbances in the derivatives of if) 
are magnified. We further observed reflections at the re- 
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fincmcnt boundaries as this pulse passes through them. 
This may explain the fluctuations in r]{x) shown in Fig. [51 
While one could continue to fine-tune a formula depen- 
dent on the conformal factor to deal with these problems, 
we looked in a different direction to determine the form 
of the damping parameter. 








FIG. 7: Form of 77(f) for a single spinning puncture sitting 
at 2; = using Eq. Q after simulation time t = 50.56 M 
(solid black line) and t = 101.25 M(dashed red line) over 
x— direction. 



B. Formulas for rj with explicit dependence on the 
position and mass of the punctures 

Since we always know the location of a puncture, and 
we know what its associated mass, we chose a form 
of damping that uses this local information throughout 
the domain. To address the demands and concerns dis- 
cussed in Section ITll and llTl A| we designed two position- 
dependent forms of r/. The two forms we tested are 



r]{r^ = A + 



Co 



l + wi(ff)" 1 + ^2(^2)"' 



and 



77(f) = A + Cie-™i (^1)" + (726-™= 



(5) 



(6) 



In Eqs. (O and ([6]), wi and W2 are required to be positive, 
unitless parameters which can be chosen to change the 
width of the functions. The power n is taken to be a 
positive integer which determines the fall-off rate. The 
constants A, Ci, and C2 are then chosen to provide the 
desired values of 77 at the punctures, and at at infinity. 



Lastly, fi and f2 are defined as 



-r| 



where i is 



either one or two, and fi is the position of the i'th black 
hole. 

The definition of fi is chosen to naturally scale the fall- 
off to the separation of the black holes, wi , W2, and n can 
be chosen to change the overall fall-off. Our work focuses 



on the choice wi = 'W2 = w and n = 1. Following [231 ] . 
we construct the damping factor to have units of inverse 
mass. We choose A = 2/Mtot, where Mtot = Mi + AI2 
is defined as the sum of the irreducible masses. We then 
take d = l/M^-A. It is then evident that both Eqs. ^ 
and ^ will give a constant value of 77 = 2/Mtot in the 
equal mass case. 

We designed the two formulas for 77 in order to test 
the value of using fundamentally different functions. In 
our simulations, we found little noticeable difference in 
the application of one compared to the other. In the 
absence of such a difference, it becomes more beneficial 
to use Eq. (O, as Gaussians are computationally more 
expensive. It should be pointed out that Eq. ([5]) is very 
similar to Eq. (13) suggested in [l^, and we believe the 
following results are very similar to what would be found 
using that form for the damping. Going into the present 
work, wc have no ansatz which might suggest these forms 
of damping yield wave forms which are any better than 
the use of any previous form of 77. However, as will be 
seen in the results sections, the waveforms we get from 
unequal mass binaries show noticeable improvement over 
the constant 77 case. 



IV. RESULTS 

For data analysis purposes, we are mainly interested in 
the properties of the emitted gravitational waves of the 
black hole binary systems under study. Hence, it is im- 
portant to check how the changes in the gauge alter the 
extracted waves. In the context of gravitational wave ex- 
traction, ^4 is only first order invariant under coordinate 
transformations. In addition, we have to chose an extrac- 
tion radius for the computation of modes, which is 
also coordinate dependent. Although the last point can 
be partly addressed by extrapolation of Tox — >■ oo, it is 
a priori not clear how much a change of coordinates af- 
fects the gravitational waves. Furthermore, a change of 
coordinates implies an effective change of the numerical 
resolution, and for practical purposes we have to ask how 
much waveforms differ at a given finite resolution. 



A. Waveform comparison using formula ([5]) 

The results in the following section refers to the use of 
Eq. ([5]). We compare numerical simulations using three 
different grid configurations, which correspond to three 
different resolutions. In the terminology of [l7| . the grid 
set-ups are 0[5x 64 : 7x128 : 5], 0[5x72 : 7x144 : 5], and 
^[5 x 80 : 7 X 160 : 5], which corresponds to resolutions on 
the finest grids of 3M/320 {N = 64), Af/120 {N = 72) 
and 3M/400 {N = 80), respectively. When referring 
to results from different resolutions, we will from here 
on use the number of grid points on the finest grid, N , 
to distinguish between them. In this subsection, we use 
wi = W2 = 12 and n ^ \ in Eq. ([S]). As test system we 
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use an unequal mass black hole binary with mass ratio 
mijm,\ = 4 and an initial separation of Z) = 5 M without 
spins in quasi-circular orbits. 

For orientation, Fig. |S] shows the amplitude of the 22- 
mode, ^22, computed with the standard gauge 77 = 2/M 
(displayed as solid lines) and with the new r](r) using 
Eq. ([5]) (displayed as non-solid lines). The three differ- 
ent colors correspond to the three resolutions. The inset 
shows a larger time range of the simulation, while the 
main plot concentrates on the time frame around merger. 
The plot gives a course view of the closeness of the results 
we obtain with standard and new gauges. 

In Fig. |9l we plot the relative differences between 
the amplitudes at low and medium (solid lines), and 
medium and high resolution (non-solid lines) obtained 
with r\ ~ 2/M (light gray lines) as well as 77 (r) (Eq. ([S])) 
(black lines). Here, we find the maximum error be- 
tween the low and medium resolution of the series using 
T] = 2/M amounts to about 12% (solid gray curve). Be- 
tween medium and high resolution (dashed gray curve), 
we find a smaller relative error, but it still goes up to 
7% at the end of the simulation. Employing Eq. ([Sj, 
the maximum amplitude error between low and medium 
resolution (solid black line) is only about 4%, and there- 
fore even smaller than the error between medium and 
high resolution for the constant damping case. Between 
medium and high resolution, the relative amplitude dif- 
ferences for Eq. ([5]) are in general smaller than the ones 
between low and medium resolution, although the maxi- 
mum error is comparable to it (dot-dashed black line). 

We repeat the previous analysis for the phase of the 
22-mode, (j)22- Again, we compare the errors between res- 
olutions in a fixed gauge. Figure fTOl shows that the error 
between lowest and medium resolution using 77 — 2/M 
(solid gray line) grows up to about 0.31 radians. For the 
differences between medium and high resolution (dashed 
line) we find a maximal error of 0.2 radians for 77 = 2/M . 
For 77(7^) following Eq. ([5]), the phase error between low 
and medium resolution is only about 0.19 radians (solid 
black line) and decreases to 0.1 radians between medium 
and high resolution (dot-dashed line). Again, employ- 
ing the position dependent form of 77, Eq. ([5]), the error 
between lowest and medium resolution is lower than the 
one we obtain for constant rj between medium and high 
resolution. The results for amplitude and phase error 
suggest that we can achieve the same accuracy with less 
computational resources using a position-dependent 77(f). 
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FIG. 8: Amplitude of the 22-mode of '^4, of a binary with 
mass ratio 4:1 and initial separation D — 5M. The different 
colors correspond to three different resolutions according to 
the grid setup described in the text. The solid lines are results 
for -q = 2/M, the dashed, dotted and dot-dashed ones are for 
r]{f) (Eq. (O). The inset shows the simulation from shortly 
after the junk radiation passed, in the main plot we zoom into 
the region of highest amplitude (near the merger) . 
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FIG. 9: Relative differences of the amplitude of the 22-mode 
of ^4 between resolutions = 64 and N — 72 (gray solid 
curve) as well as N = 72 and N = SO (gray dashed curve) 
when using rj — 2/M. The same for rj{f) (Eq. (O) between 
A = 64 and N = 72 (black solid curve) and N = 72 and 
N = SQ (black dot-dashed curve). The physical situation is 
the same as in Fig. [8l The maximum differences are above 
10%, comparing low and medium resolution of the constant 
?7 simulations (gray solid line). 



B. Waveform comparison using formula (f6)) 

We repeated the analysis of Sec. IIV Al with the wave- 
forms we obtain using Eq. ([6]) (with wi — W2 — ^2 and 
n = 1). We use the same initial conditions (mass ratio 
4 : 1, D = 5M, no spins), and compare the amplitudes 
and phases of the 22-mode of ^4 with the results of the 
77 = 2/M-runs. The grid configurations remain the same. 

The results are very similar to the ones we obtained in 



Figs. [9I and [TOl and we therefore do not show them here. 
Although Eqs. ([5]) and © result in different shapes for 
77(7^), is very similar. Therefore, the comparison to 
77 = 2/M naturally gives very similar results, too. The 
phase differences between results from Eqs. ([5|) and ([6]) 
at a given resolution are shown in Fig. [TT] These are, 
with a maximum phase error of 0.004 radians, very small 
compared to the phase errors between resolutions, which, 
at minimum, are about 0.1 radian (see Fig. [T0\\ . Fig. [T2l 
compares the phase error between low and medium (solid 
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FIG. 10: Phase difTerences between lowest and medium reso- 
lution for the series using rj — 2/M (solid gray line) and ri{f) 
(Eq. ((Sjl) (solid black line) as well as between medium and 
high resolution for ry — 2/M (dashed gray line) and for rj{f) 
(Eq. ([5])) (dot-dashed black line). The physical situation is 
the one of Fig. [S] 
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FIG. 11: Phase differences between waveforms obtained 
with Eq. ((Sjl and Eq. ([6]) in three different resolutions (solid, 
dashed, dotted-dashed lines) for mass ratio 4:1, D — 5M. 



lines), and medium and high resolution (dotted-dashed 
and dashed line) of Eq. ([6]) (gray) to the ones of Eq. ([5]) 
(black). For comparison, the error between medium and 
high resolution is also plotted for Eq. (|4]) in this figure 
(dotted line). The plot indicates that the errors between 
resolutions are in good agreement for the different posi- 
tion dependent formulas of rj. 



FIG. 12: Phase difference between waveforms at low and 
medium resolution (solid lines) and medium and high reso- 
lution (dotted-dashed and dashed line) using either Eq. ([5]) 
(black lines) or Eq. ^ (gray lines) for mass ratio 4:1, D = 
5 M. For comparison, we also show the phase difference ob- 
tained with Eq. (|4]) between medium and high resolution (dot- 
ted line). 
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FIG. 13: i-component of the shift vector in a--direction after 
160 M of evolution of the system with mass ratio 4 : 1 and 
D — 10 M. The black, dot-dashed line refers to the use of a 
constant damping rj, while the black, solid line uses Eq. 
The gray, dashed line is for the use of Eq. ((6| and the gray, 
dotted one for Eq. ((5|. Except for the constant rj (black, dot- 
dashed line), the results in this plot are indistinguishable. 



C. Behavior of the shift vector 

In [23I I , we found an unusual behavior of the shift vec- 
tor. This is illustrated in Fig. [T31 where we plot the 
x-component of the shift, /?^, in the cc-direction after 
160 M of evolution (this means approximately 80 M af- 
ter merger) for all four versions of the damping constant 
we used for comparison in this paper before, and for the 
same binary configuration as the one used in Sees. IIV Al 



and If V B] Like in [2^, we find that using Eq. (U) results 
in a shift which falls off to zero too slowly towards the 
outer boundary, and which develops a "bump" (black, 
solid line), while the constant damping case (black, dot- 
dashed line) falls off to zero quickly. Employing Eqs. ([5]) 
or ^ avoids this undesirable feature. After merger, the 
shift falls off to zero when going away from the punctures 
as it does in the constant damping case (gray dashed and 
dotted lines). Using Eq. ©or ((51) prevents unwanted co- 
ordinate drifts at the end of the simulations. 
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FIG. 14: Shown is the time dependence of the ratio between 
the coordinate areas of the apparent horizons of both black 
holes in a simulation with mass ratio 4 : 1 with initial separa- 
tion D = 5M. The black, blue and red lines use ri{r), Eq. ^ 
with varying values of the width parameter w. The orange hne 
(dash-dot-dot) uses the constant damping 77 = 2/A'I and the 
green (dash-dot) one refers to the result of [23| with Eq. Q. 
Using Eq. ((SJ, the coordinate areas can be varied with re- 
spect to each other depending on the choice of w. A ratio of 
1 means the black holes have the same size on the numerical 
grid. 



V. DISCUSSION 



the numerical grid. Fig. [14] plots the ratio of the grid- 
area of larger apparent horizon to the smaller apparent 
horizon as a function of time for w- values of 0.1, 0.5 and 
for 200, all with n = 1. Also plotted is the relative co- 
ordinate size for the same binaries using a constant 77 
in dashed, double-dotted line, and for using Eq. ^ in 
a blue dashcd-dotted line. All the evolutions show an 
immediate dip, and then increase in the grid-area ratio 
during the course of the evolution. While a very low ra- 
tio was found using Eq. the orange dotted line was 
later found for the choices of n = 3 with wi = 0.01 and 
W2 = 0.0001 with Eq. ([5]). Due to this freedom in the 
implementation of our explicit formula for the damping, 
it may be possible to further reduce the relative grid size 
of the black holes. This effect could be important in eas- 
ing the computational difficulty of running a numerical 
simulation for unequal mass binaries. 

Having a form of 77 that leads to stable evolutions for 
any mass-ratio is an important step towards the numer- 
ical evolution of binary black holes in the intermediate 
mass-ratio. We believe the form given in Eq. ([5]) pro- 
vides such a damping factor at a low computational cost, 
although the test results presented are limited to mass 
ratio 4:1. We plan to examine larger mass ratios in 
future work. The new method should allow binary sim- 
ulations for mass ratio 10 : 1, or even 100 : 1. It remains 
to be seen whether other issues than the gauge are now 
the limiting factor for simulations at large mass ratios. 



In this work, we examined the role that the damping 
factor, r], plays in the evolution of the shift when using 
the gamma driver. In particular, we examined the range 
of values allowed in various evolutions, and what effects 
showed up because of the value chosen. We then de- 
signed a form of ry for the evolution of binary black holes 
which provides appropriate values both near the individ- 
ual punctures and far away from them with a smooth 
transition in between. 

In Sec. lIVl we directly examined the waveforms for the 
case using Eq. ([S]), where wi = W2 = 12 and n = 1. While 
the form of 77 is predictable, and can be easily adjusted 
for stability, we also saw that the waveforms produced 
using this definition showed less deviation with increas- 
ing resolution than using a constant 77. When examining 
the waveforms produced using Eq. we found simi- 
lar results. In the absence of a noticeable difference in 
the quality of the waveforms, Eq. ([5|) is computationally 
cheaper, and, as such, is our preferred definition for the 
damping. 

We have already pointed out a certain freedom to 
pick parameters in Eqs. ^ and We did perform 

some experimentation along this line where we varied 
w ^ wi ^ W2 to see if we could get a useful effect of the 
coordinate size of the apparent horizons on the numer- 
ical grid. In [13, H^, it was noticed that the damping 
coefficient affects the coordinate location of the apparent 
horizon, and therefore the resolution of the black hole on 
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